/**
    \file  kmeans.cc
    \brief Driver program to test and profile libgist's k-means.

    This program test-drives libgist's implementation of the k-means
    clustering algorithm. It can be used to collect performance metrics
    and also compare libgist's k-means implementation to OpenCV (and,
    possibly, other) k-means implementations.

    NOTE: Although it is located in the test subdirectory of the libgist
    codebase, this program is not actually a unit test to ascertain
    correctness. Rather, it is a test to gauge clustering performance.

*/

/*
   This file is part of libgist.

   libgist is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2 of the License, or (at your
   option) any later version.

   libgist is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public License
   along with libgist; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
*/

/*
   REVISION HISTORY

   $HeadURL: http://libgist.googlecode.com/svn/branches/lazebnik_dev/test/kmeans.cc $
   $Id: kmeans.cc 179 2012-04-12 03:18:05Z libgist@gmail.com $
*/

//------------------------------ HEADERS --------------------------------

// libgist
#include "kmeans.hh"

// Boost
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>

#include <boost/random/variate_generator.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/mersenne_twister.hpp>

#include <boost/date_time/posix_time/posix_time_types.hpp>
#include <boost/program_options.hpp>
#include <boost/lambda/lambda.hpp>

// Standard C++
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <vector>
#include <iterator>
#include <stdexcept>
#include <string>
#include <utility>
#include <limits>

// Standard C
#include <time.h>
#include <stdlib.h>

//----------------------- PROGRAM RETURN CODES --------------------------

namespace {

enum error_codes {
   BAD_OPTION = 1,
} ;

} // end of local anonymous namespace encapsulating above constants

//---------------------------- BASIC HELP -------------------------------

namespace {

std::string purpose()
{
   return "kmeans -- test-drive libgist's implementation of k-means" ;
}

std::string copyright()
{
   return "Copyright (C) 2012 Manu Viswanathan" ;
}

std::string usage()
{
   const char* usage =
"usage: kmeans -g [-n  points ] [-d dimensions] [-r min max ] data_file\n"
"       kmeans -o [-n  points ] [-d dimensions] [-k clusters] \\\n"
"                 [-x max_iter] data_file\n"
"       kmeans    [-n  points ] [-d dimensions] [-k clusters] \\\n"
"                 [-x max_iter] [-t threads] data_file\n"
"       kmeans [-h|--help]\n"
"       kmeans [-w|--warranty]\n" ;

   return usage ;
}

void warranty()
{
   std::cout << purpose()   << '\n' ;
   std::cout << copyright() << "\n\n" ;

   const char* gpl_notice =
"libgist is free software; you can redistribute it and/or modify it under the\n"
"terms of the GNU General Public License as published by the Free Software\n"
"Foundation; either version 2 of the License, or (at your option) any\n"
"later version.\n\n"

"libgist is distributed in the hope that it will be useful, but WITHOUT ANY\n"
"WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS\n"
"FOR A PARTICULAR PURPOSE. See the GNU General Public License for more\n"
"details.\n\n"

"You should have received a copy of the GNU General Public License along\n"
"with libgist; if not, write to the Free Software Foundation, Inc., 59 Temple\n"
"Place, Suite 330, Boston, MA 02111-1307 USA.\n" ;
   std::cout << gpl_notice ;
}

void help()
{
   std::cout << purpose()  << '\n' ;
   std::cout << copyright()<< "\n\n" ;
   std::cout << usage()    << '\n' ;

   const char* help =
"This program test-drives libgist's implementation of the k-means clustering\n"
"algorithm. It is meant to be used to collect performance data and compare\n"
"with OpenCV's implementation of k-means.\n\n"

"This program operates in one of two modes, viz., either data generation or\n"
"clustering. In data generation mode, it will generate random data points of\n"
"the specified dimensionality and save this data to the specified file. In\n"
"clustering mode, it will read data from a previously generated file and\n"
"cluster the data using either libgist's native, multithreaded k-means or\n"
"OpenCV's k-means.\n\n"

"Here is a brief explanation of the supported command line options:\n\n"

"-g|--generate-data\n"
"\tThis flag turns on data generation mode. Without it, the program\n"
"\tworks in clustering mode.\n\n"

"-r|--range min max\n"
"\tThis option can be used in data generation mode to specify the\n"
"\trange of data values for each dimension. The default range is\n"
"\t[0, 128].\n\n"

"-n|--points n\n"
"\tThis option specifies the number of data points. In data generation\n"
"\tmode, it will be the number of data points to generate. In clustering\n"
"\tmode, it specifies the number of data points to read from the input\n"
"\tfile.\n\n"

"\tThis is an optional parameter. If not specified, it defaults to a\n"
"\tvalue of one million. Thus, in data generation mode, by default,\n"
"\tthe program will generate one million points. And in clustering\n"
"\tmode, it will (try to) read one million points from the data file.\n\n"

"-d|--dimensions d\n"
"\tThis option specifies the number of dimensions in the input data.\n"
"\tIt defaults to a value of 128; i.e., if it is not specified, the\n"
"\tprogram will assume it is working with 128-dimensional data.\n\n"

"-k|--clusters k\n"
"\tThis option specifies the number of clusters to create. It is only\n"
"\tvalid in clustering mode and its default value is 200.\n\n"

"-x|--max-iterations x\n"
"\tThis option specifies the maximum number of iterations after which\n"
"\tclustering should terminate. If the algorithm converges before\n"
"\tgoing through these many iterations, it will terminate earlier.\n"
"\tThe default value for this parameter is 5000.\n\n"

"-t|--threads t\n"
"\tlibgist's implementation of k-means is multithreaded. This option\n"
"\tcan be used to specify the number of threads to use while\n"
"\tclustering with libgist's native k-means. If not specified, it\n"
"\tdefaults to one, i.e., clustering will be single-threaded.\n\n"

"-o|--opencv\n"
"\tIn clustering mode, by default, this program will use libgist's\n"
"\timplementation of k-means. However, with this flag, it will use\n"
"\tOpenCV's k-means instead. Note that you cannot use -t when\n"
"\tclustering with OpenCV.\n\n"

"-h|--help\n"
"\tPrint this (hopefully) useful info.\n\n"

"-w|--warranty\n"
"\tPrint warranty (or rather, lack thereof).\n\n" ;
   std::cout << help ;
}

} // end of local anonymous namespace encapsulating above helpers

//----------------------- COMMAND LINE PARSING --------------------------

namespace {

// Quick helper to package command line args together
struct cmdline {
   bool generate_data, opencv ;
   int  n, d, k, x, t ;
   std::pair<float, float> range ;
   std::string data_file;

   cmdline() ;
} ;

// Init command line args to their default values
cmdline::cmdline()
   : generate_data(false), opencv(false),
     n(1000000), d(128), k(200), x(5000), t(1),
     range(0.0f, 128.0f)
{}

// Helper function to report an error message related to command line
// options and exit program with appropriate error code.
void option_error(const std::string& msg)
{
   std::cerr << "kmeans: " << msg << '\n' ;
   std::cerr << usage() ;
   exit(BAD_OPTION) ;
}

// Helper function to take care of the annoying details of using
// Boost.program_options to get at the command line arguments.
//
// DEVNOTE: This code is lifted almost verbatim from the
// Boost.program_options tutorial and examples that comes with the Boost
// documentation. There may be better (i.e., neater, more effective,
// clearer, more efficient, whatever) ways to use the library.
cmdline parse(int argc, char* argv[])
{
   namespace po = boost::program_options ;
   cmdline params ;
   try
   {
      std::vector<float> range(2) ;
      range[0] = std::numeric_limits<float>::min() ;
      range[1] = std::numeric_limits<float>::max() ;

      // Specify the command line options
      po::options_description options("Command line options") ;
      options.add_options()
         ("generate-data,g",
          po::bool_switch(&params.generate_data)->default_value(false),
          "turn on data generation mode")

         ("range,r", po::value<std::vector<float> >(&range)->multitoken(),
          "range of values for each data dimension")

         ("points,n", po::value<int>(&params.n), "number of data points")

         ("dimensions,d", po::value<int>(&params.d),
          "number of dimensions in each data point")

         ("clusters,k", po::value<int>(&params.k),
          "number of clusters to create")

         ("max-iterations,x", po::value<int>(&params.x),
          "maximum number of iterations after which clustering terminates")

         ("threads,t", po::value<int>(&params.t),
          "number of threads to use when clustering with libgist's "
          "native, multithreaded k-means")

         ("opencv,o", po::bool_switch(&params.opencv)->default_value(false),
          "cluster with OpenCV rather than libgist")

         ("data-file,f", po::value<std::string>(&params.data_file),
          "name of the data file")

         ("help,h", "print some help")
         ("warranty,w", "print warranty") ;

      // Convert first positional command line argument to -f
      po::positional_options_description p ;
      p.add("data-file", 1) ;

      // Setup done: now parse argc and argv...
      po::variables_map vm ;
      po::store(po::command_line_parser(argc, argv).
                options(options).positional(p).run(), vm) ;
      po::notify(vm) ;

      // Process -h and -w
      if (vm.count("help")) {
         help() ;
         exit(0);
      }
      if (vm.count("warranty")) {
         warranty() ;
         exit(0) ;
      }

      // Mode-based sanity checks
      if (params.generate_data) // data generation mode
      {
         if (vm.count("clusters")       || vm.count("threads") ||
             vm.count("max-iterations") || params.opencv)
            option_error("cannot use -k, -o, -x, or -t with -g") ;

         if (vm.count("range")) {
            if (range[0] > range[1])
               std::swap(range[0], range[1]) ;
            params.range = std::make_pair(range[0], range[1]) ;
         }
      }
      else // clustering mode
      {
         if (vm.count("range"))
            option_error("cannot use -r in clustering mode") ;

         if (params.opencv && vm.count("threads"))
            option_error("cannot use -t with -o") ;
      }

      // Sanity checks common to both modes
      if ((vm.count("points")         && params.n < 1) ||
          (vm.count("dimensions")     && params.d < 1) ||
          (vm.count("clusters")       && params.k < 1) ||
          (vm.count("max-iterations") && params.x < 1) ||
          (vm.count("threads")        && params.t < 1))
         option_error("-n, -d, -k, -x, and -t require positive numbers") ;

      if (params.data_file.empty())
         option_error("please specify data file name") ;
   }
   catch (po::error& e)
   {
      option_error(e.what()) ;
   }
   return params ;
}

} // end of local anonymous namespace encapsulating above helpers

//----------------------- I/O AND DEBUG SUPPORT -------------------------

namespace {

// Quick helper to dump STL vectors to stdout
template<typename T>
void dump(const std::vector<T>& v)
{
   std::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " ")) ;
   std::cout << '\n' ;
}

// Quick helper to dump Boost.uBLAS matrices to stdout
template<typename T>
void dump(const boost::numeric::ublas::matrix<T>& m, const std::string& prefix)
{
   const int d = m.size1() ;
   const int n = m.size2() ;
   std::cout << prefix << ": " << d << 'x' << n << '\n' ;
   for (int i = 0; i < d; ++i) {
      boost::numeric::ublas::vector<T> v = boost::numeric::ublas::row(m, i) ;
      std::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " ")) ;
      std::cout << '\n' ;
   }
}

// Quick helper to read back the data file to ensure that serialization
// worked properly.
void dump(int n, int d, const std::string& file_name)
{
   std::cout << "reading " << n << 'x' << d << '\n' ;
   std::ifstream file(file_name.c_str(), std::ios::binary) ;
   const int D = d * sizeof(float) ;
   std::vector<float> numbers(d) ;
   for (int i = 0; i < n; ++i) {
      file.read(reinterpret_cast<char*>(&numbers[0]), D) ;
      dump(numbers) ;
   }
}

// Read previously generated data file and return a d-by-n matrix
template<typename T>
boost::numeric::ublas::matrix<T>
read(int d, int n, const std::string& file_name)
{
   boost::numeric::ublas::matrix<T> m(d, n) ;
   std::ifstream file(file_name.c_str(), std::ios::binary) ;
   const int D = d * sizeof(T) ;
   std::vector<T> v(d) ;
   for (int i = 0; i < n; ++i) {
      file.read(reinterpret_cast<char*>(&v[0]), D) ;
      if (!file.good()) {
         std::ostringstream msg ;
         msg << file_name << ": not enough data for "
             << d << 'x'  << n << " matrix" ;
         throw std::runtime_error(msg.str());
      }
      std::copy(v.begin(), v.end(),
                boost::numeric::ublas::column(m, i).begin()) ;
   }
   return m ;
}

} // end of local anonymous namespace encapsulating above helpers

//-------------------------- DATA GENERATION ----------------------------

namespace {

// Generate n random d-dimensional points with each dimension's value
// lying in the range [min, max) and store the results to the specified
// file.
//
// NOTE: The data file does not encode/impose any structure on the
// numbers it generates. That is, the data file is not a d-by-n matrix or
// any such thing; rather it is a simple flat sequence of d*n numbers
// stored in binary form. Functions that read this file can impose
// whatever structure they want.
void
generate_data(int n, int d, float min, float max, const std::string& file_name)
{
   typedef boost::mt19937 gen_t ;
   typedef boost::uniform_real<float> dist_t ;

   gen_t  gen(static_cast<unsigned int>(time(0))) ;
   dist_t dist(min, max) ;
   boost::variate_generator<gen_t&, dist_t> rng(gen, dist) ;

   std::ofstream file(file_name.c_str(), std::ios::binary) ;
   //std::cout << "writing " << n << 'x' << d << '\n' ;

   const int D = d * sizeof(float) ;
   std::vector<float> numbers(d) ;
   for (int i = 0; i < n; ++i) {
      std::generate(numbers.begin(), numbers.end(), rng) ;
      file.write(reinterpret_cast<const char*>(&numbers[0]), D) ;
      //dump(numbers) ;
   }
}

} // end of local anonymous namespace encapsulating above helpers

//---------------------- CLUSTERING WITH LIBGIST ------------------------

namespace {

void kmeans(int n, int d, int k, int x, int t, const std::string& file_name)
{
   typedef gist::kmeans<float> km ;
   try
   {
      boost::numeric::ublas::matrix<float> data(read<float>(d, n, file_name)) ;
      //std::cout << "want to read " << d << 'x' << n << '\n' ;
      //dump(data, "read") ;

      using boost::posix_time::ptime ;
      using boost::posix_time::microsec_clock ;
      ptime start(microsec_clock::local_time()) ;

      km clustering(t, x) ;
      clustering.centroids(k, data) ;
      ptime finis(microsec_clock::local_time()) ;
      float time_taken = (finis - start).total_milliseconds()/1000.0f ;

      std::cout << n << ' ' << k << ' ' << d << ' ' << t << ' '
                << clustering.iterations()   << ' ' << time_taken << ' '
                << clustering.wss() << '\n' ;
   }
   catch (gist::kmeans_errors::empty_cluster& e)
   {
      std::cerr << "kmeans: empty cluster " << e.cluster
                << " at iteration " << e.iteration << " when clustering "
                << n << ' ' << d << "-dimensional points into "
                << k << " clusters\n" ;
   }
}

} // end of local anonymous namespace encapsulating above helpers

//------------------------------- MAIN ----------------------------------

int main(int argc, char* argv[])
{
   try
   {
      cmdline params = parse(argc, argv) ;
      if (params.generate_data) {
         generate_data(params.n, params.d,
                       params.range.first, params.range.second,
                       params.data_file) ;
         //dump(params.n, params.d, params.data_file) ;
      }
      else if (params.opencv)
         std::cout << "num points: "     << params.n << '\n'
                   << "num clusters: "   << params.k << '\n'
                   << "num dimensions: " << params.d << '\n'
                   << "max iterations: " << params.x << '\n'
                   << "num threads: "    << params.t << '\n'
                   << "clustering library: "
                   << (params.opencv ? "OpenCV" : "libgist") << '\n' ;
      else // libgist kmeans ==> can be multithreaded
         kmeans(params.n, params.d,
                params.k, params.x, params.t, params.data_file) ;
   }
   catch (std::exception& e)
   {
      std::cerr << "kmeans: " << e.what() << '\n' ;
      return 255 ;
   }
   return 0 ;
}

//-----------------------------------------------------------------------

/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */
